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Abstract 

Prognostics is centered on predicting the time of and time un- 
til adverse events in components, subsystems, and systems. 
It typically involves both a state estimation phase, in which 
the current health state of a system is identified, and a pre- 
diction phase, in which the state is projected forward in time. 
Since prognostics is mainly a prediction problem, prognos- 
tic approaches cannot avoid uncertainty, which arises due to 
several sources. Prognostics algorithms must both character- 
ize this uncertainty and incorporate it into the predictions so 
that informed decisions can be made about the system. In this 
paper, we describe three methods to solve these problems, in- 
cluding Monte Carlo-, unscented transform-, and first-order 
reliability-based methods. Using a planetary rover as a case 
study, we demonstrate and compare the different methods in 
simulation for battery end-of-discharge prediction. 

1. Introduction 

Prognostics focuses on predicting the time of and time un- 
til adverse events in components, subsystems, and systems. 
Model-based methods consist of an estimation phase, in 
which the current health state of a system is identified, fol- 
lowed by a prediction phase, in which the state is simulated 
until the adverse event. Because prognostics is mainly a pre- 
diction problem, uncertainty, which manifests due to many 
sources, cannot be avoided. This uncertainty will determine 
how the system evolves in time, i.e., the system evolution is 
a random process. To approach this problem in a system- 
atic way, there are two problems to solve: (i) characteriz- 
ing the sources of uncertainty, and ( ii ) quantifying the com- 
bined effect of the different sources of uncertainty within the 
prediction. With proper estimation of the prediction uncer- 
tainty, predictions can then be used to make informed deci- 
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sions about the system. 

In many applications, the largest source of uncertainty is that 
associated with the future inputs, and so most work in uncer- 
tainty quantification for prognostics has focused in that area. 
In the context of fatigue damage prognosis, different types of 
methods (Ling et al., 2011) such as rainflow counting tech- 
niques, auto-regressive moving-average models, and Markov 
processes have been used for constructing future loading tra- 
jectories. In (Sankararaman, Ling, Shantz, & Mahadevan, 
2011), the authors construct future input trajectories as se- 
quential blocks of constant-amplitude loading, and such tra- 
jectories are sampled in the prediction algorithm. In (Saha, 
Quach, & Goebel, 2012) the authors characterize the future 
inputs by determining statistics of the battery loading for typ- 
ical unmanned aerial vehicle maneuvers based on past flight 
data, and construct future input trajectories as constrained se- 
quences of flight maneuvers. Constant loading is assumed 
in (Luo, Pattipati, Qiao, & Chigusa, 2008) for a vehicle sus- 
pension system, and predictions are made for a weighted set 
of three different loading values. 

Once each source of uncertainty has been characterized, it 
must be accounted for within the prediction, and thereby their 
combined effect on the overall prediction must be quanti- 
fied. In previous work (Daigle, Saxena, & Goebel, 2012), 
we investigated methods for accouting for future input uncer- 
tainty in the predictions and introduced the unscented trans- 
form (UT) method for efficiently estimating the future input 
uncertainty, however, methods for future input characteriza- 
tion were not discussed and only constant-amplitude loading 
was considered. In other work (Sankararaman, Daigle, Sax- 
ena, & Goebel, 2013; Sankararaman & Goebel, 2013), we 
investigated the use of analytical methods, namely, first-order 
reliability (FORM) based methods for propagating the future 
input uncertainty, however it also was limited to constant- 
amplitude loading, and the methods were applied only for 
offline analysis and not within online prognostic algorithms. 
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In this paper, we adopt a model-based prognostics approach 
(Orchard & Vachtsevanos, 2009; Daigle & Goebel, 2013; 
Luo et al., 2008). We discuss the uncertainty representa- 
tion and quantification problem, and develop a novel gen- 
eralized framework using the notion of surrogate variables, 
allowing the representation of state uncertainty, parame- 
ter uncertainty, future input uncertainty, and process noise 
in a common framework that allows consideration of both 
constant- and variable -amplitude loading situations. We ex- 
plore three methods for prediction with uncertainty quantifi- 
cation, namely, Monte Carlo sampling, UT sampling, and in- 
verse FORM. Using a rover battery system as a case study, we 
describe two methods for future input uncertainty characteri- 
zation, implement the prediction methods, and compare their 
performance for battery end-of-discharge prediction in simu- 
lated constant- and variable-amplitude loading scenarios. 

The paper is organized as follows. In Section 2, we define the 
prognostics problem and develop the uncertainty representa- 
tion framework. In Section 3, we describe methods for han- 
dling uncertainty within the prediction process. In Section 4, 
we introduce the rover case study and present several exam- 
ples in simulation demonstrating the methods and comparing 
their performance. Section 5 concludes the paper. 

2. Model-based Prognostics 

We first formulate the prognostics problem, and develop the 
uncertainty representation framework. We then provide an 
architecture for model-based prognostics. 

2.1. Problem Formulation 

We assume the system model may be generally defined as 

x(fc + 1) = f (k, x(fc), 9(k), u (fc), v(fc)), (1) 

yO) = h(fc, x(fc), 6[k), u(fc), n(fc)), (2) 

where k is the discrete time variable, x(fc) £ M” 1 ' is the 
state vector, 9(k) £ is the unknown parameter vector, 
u(fc) £ K" 1 * is the input vector, v(fc) £ K ra ” is the process 
noise vector, f is the state equation, y(k) £ K™ 1 ' is the output 

vector, n(fc) £ is the measurement noise vector, and h 

is the output equation. 1 The unknown parameter vector 9(k) 
is used to capture explicit model parameters whose values are 
unknown and time-varying stochastically. 

In prognostics, we are interested in predicting the occurrence 
of some (desirable or undesirable) event E that is defined with 
respect to the states, parameters, and inputs of the system. 
We define the event as the earliest instant that some event 
threshold T E : K n * x R" 8 x -A B, where B = {0, 1} 
changes from the value 0 to 1. That is, the time of the event 


1 Here, we use bold typeface to denote vectors, and use n a to denote the 

length of a vector a. 


kp at some time of prediction kp is defined as 

k E {k P ) = inf{fc G N; k > kp A T E (pc(k), 9{k), u(fc)) = 1}. 

(3) 

The time remaining until that event, A k E , is defined as 

Ak E (kp) A k E (kp) — kp. (4) 

In the context of systems health management, the event E 
corresponds to some undesirable event such as the failure of 
a system, some process variable out-of-range, or a similar 
type of situation. T E is defined via a set of performance con- 
straints that define what the acceptable states of the system 
are, based on x(fc), 9(k), and u (k) (Daigle & Goebel, 2013). 
In this case, k E is then conventionally called end of life, while 
A k E is conventionally called remaining useful life. 

The system evolution is a random process due to the pro- 
cess noise v(fc) and nondeterministic inputs u (k). Since the 
system evolution is a random process, at time kp. the sys- 
tem takes some trajectory out of many possible trajectories, 
therefore, k E and A k E are random variables. So, the prog- 
nostics problem is to predict the probability distribution of 
k E (Daigle, Saxena, & Goebel, 2012; Sankararaman et al., 
2013). 

Problem (Prognostics). The prognostics problem is to, at 
prediction time kp, compute p(k E (kp)\y(ko-kp)) and/or 
p(Ak E (kp)\y(k 0 :k P )). 

In practice, obtaining this distribution is difficult because the 
state at kp, system model, process noise distribution, and fu- 
ture input distribution may not be known precisely. 

2.2. Representing Uncertainty 

In order to predict k E , four sources of uncertainty must be 
dealt with: (i) the initial state at time kp, x(fc p ); (ii) the pa- 
rameter values 9(k) for all k > kp, denoted as &k P (the 
subscript kp indicates the start time of the trajectory); (Hi) 
the inputs u(k) for all k > kp, denoted as Ufc p ; and (iv) the 
process noise v(fc) for all k > kp, denoted as V/. (> . In or- 
der to make a prediction that accounts for this uncertainty, we 
require the probability distributions p(x), p(&k P ), p(Uk P ), 
and p(Vfcp). 

For describing the probability distribution of a generic trajec- 
tory A f., we introduce a set of surrogate random variables 
\ a = [A*A^ . . .]. We describe a trajectory using \ a and in- 
stead define p(A a ), which in turn defines p(A/, ). These sur- 
rogate variables can be used to describe trajectories in a vari- 
ety of ways. For example, we can associate A k(k) with A*, 
A k(k + 1) with Aq, etc. Or, we can describe a trajectory with 
a parameterized function, where the parameters are the ran- 
dom variables, e.g., A k(k) = A* + A^ sinfc. The use of the 
surrogate variables provides flexibility to the user in defining 
the distribution of trajectories. So, for the parameter, input. 
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Figure 1. Prognostics architecture. 


and process noise trajectories we introduce the surrogate vari- 
ables A g, A„, and A,,. 

In the model-based prognostics paradigm, the probability 
distribution for the initial state at the time of prediction 
kp , p(x(kp),6(kp)), is specified by an estimator, such as 
the Kalman filter, unscented Kalman filter (UKF) (Julier & 
Uhlmann, 1997), or particle filter (Arulampalam, Maskell, 
Gordon, & Clapp, 2002). The distribution maybe represented 
by a set of statistical moments like a mean vector and covari- 
ance matrix (as in the Kalman filter), or a set of weighted 
samples (as in the UKF and particle filter). This defines only 
the parameter vector at kp and not for future time steps. For 
the future values of 9 , these are drawn from a distribution 
specified by A g. In the simplest case, Aj can define 9{kp ), 
Ag can define 9{kp + 1), etc, where the distribution of each 
A l g is the same and specified by p(9{k p )) as determined by the 
estimator. Alternatively, it can be assumed that 0(k) is con- 
stant, in which case only one surrogate variable is required, 
with the distribution specified by the estimator. 

For process noise, we define p(V k ) using A„. It is often as- 
sumed in practice that, at each time instant, there is a single 
probability distribution from which the process noise values 
are drawn. The distribution is defined a priori based on the 
understanding of the system and its model. In the simplest 
case, there is a random variable for every time step, i.e. A, 1 , 
defines v(fcp), Xf, defines v{kp + 1), etc. Since the num- 
ber of random variables depends on the number of time steps, 
such a treatment potentially leads to dimensionality issues. 
Sankararaman and Goebel (Sankararaman & Goebel, 2013) 
have demonstrated that it is possible to approximate the effect 
of this process noise using an equivalent time-invariant pro- 
cess noise, i.e., a single random variable that defines a con- 
stant value of process noise for all k. In this case. A,, would 
contain only that single random variable, whose probability 
distribution will be constructed by matching the likelihood of 
occurrence of the original time-varying process noise. 

For the future input trajectories, the distribution depends on 
the particular system being prognosed and the environment 
it is operating within. As with the other trajectories, we de- 
scribe p(Uk P ) using A,,. Often, there is some knowledge 
about what the future input will be and only a few random 
variables are needed in A„. For example, in a constant- 


loading scenario, the inputs can be defined with u (k) = A* 
for k > kp. Any arbitrary function parameterized by a 
set of random variables may be used to define u (k), e.g., 
u (k) = Xi ■ sin(k), or u (k) = A* + A^ • k. The variables in 
A 0 may or may not be independent. 

To predict Ate at time kp, we require the initial state, 
x(fcp); the parameter trajectory up to kp, & kp = 
[9(kp), . . . , 9{kp)\, the process noise trajectory up to kp, 
V kp = [v (kp v(fcp)]; and the input trajectory up to 
kp, U kp = [u(fcp), . . . , u(Ap)]. With this information, kp 
can be computed simply by simulating the state forward in 
time, using the state equation f, until Tp evaluates to 1, at 
which point the current time is kp. Because we have only 
probability distributions, we need to propagate the uncer- 
tainty through this procedure in order to obtain the distri- 
bution for kp (Sankararaman et ah, 2013; Sankararaman & 
Goebel, 2013). Methods to do such uncertainty propagation 
will be described in Section 3. 

2.3. Prognostics Architecture 

We adopt a model-based prognostics architecture (Daigle & 
Goebel, 2013), in which there are two sequential problems, 
(/) the estimation problem, which requires determining a 
joint state-parameter estimate p(x(fc), 9(k)\y(ko-kp)) based 
on the history of observations up to time k, y(ko:kp), and 
(n) the prediction problem, which determines at kp, using 
p{x(k), 9(k)\y(k 0 :k P )), p( Xg), p(X u ), and p( A„), a proba- 
bility distribution p(kp(kp)\y(ko:kp)). The distribution for 
Akp can be trivially computed fromp(kp(kp)ly(ko:kp)) by 
subtracting kp from kp(kp). 

The prognostics architecture is shown in Fig. 1. In discrete 
time k, the system is provided with inputs u k and provides 
measured outputs y k . The estimation module uses this in- 
formation, along with the system model, to compute an es- 
timate p(x(fc), 9(k)\y(ko'.k)). The prediction module uses 
the joint state-parameter distribution and the system model, 
along with the distributions for the surrogate variables, p(\g), 
p(X u ), and p(X v ), to compute the probability distribution 
p{kp(kp)\y(ko'.kp)) at given prediction times kp. 
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Algorithm 1 k E (k P ) <- P(x(fc P ), Q kp , U fcp , V fcp ) 

1 1 k t — k p 
2: x(fc) -t— x(fcp) 

3: while T E (x(k),& kp (k),\Jk P (k)) = 0 do 
4: x(fc + 1) <- f (fc, x(fc), © fep (fc), U fcp (fc), V fep (fc)) 

5: k i — k -{- 1 

6 : x(fc) <— x(fc + 1) 

7: end while 

8 : kE(kp) < — k 


3. Prediction 

Prediction is initiated at a given time kp using the current 
joint state-parameter estimate, p(x(kp), 9(kp)\y(ko'.kp)). 
Approaches to determine this estimate are reviewed 
in (Daigle, Saha, & Goebel, 2012) and are not described here. 
The goal is to compute p(kp (fcp)|y (ko'-kp)) using the state- 
parameter estimates and assumptions about uncertainty re- 
garding the future parameter, input, and process noise values. 

Consider one realization of each of the uncertain quantities at 
prediction time kp: the state x(fc P ), the parameter trajectory 
&k P , the input trajectory U; Cp , and the process noise trajec- 
tory V k p . Then, the corresponding realization of kp can be 
computed with the system model as shown in Algorithm 1 . In 
Algorithm 1, the function P simulates the system model until 
the threshold Tp evaluates to 1. 

This algorithm requires computing first realizations of the 
state-parameter distribution, the parameter trajectory, the in- 
put trajectory, and the process noise trajectory. As described 
in Section 2, the distribution for the state comes from an es- 
timator, and the distributions for the parameter, input, and 
process noise trajectories are defined indirectly by a set of 
siurogate variables. At the higher level, we are interested in 
computing the distribution for kp from the distributions for 
p(x(fc P ), 9(k P )), p(Xg), p( X u ), and p(X v ). 

The function that takes these siurogate variables and com- 
putes a distribution for kp, which we refer to as G, 
is the real function we are interested in, i.e., p(kp) = 
G(p(x(k P )),p(9(k P )),p(X s ),p(X u ),p(X v )). To compute 
p(kp), we must propagate the uncertainty through this func- 
tion. That is, predicting p(kp(kp)\y(ko:kp)) is an uncer- 
tainty propagation problem. In the following subsections, 
we describe three different methods with which to solve this 
problem. They each compute realizations of the state, param- 
eter trajectory, input trajectory, and process noise trajectory, 
and call the P function to obtain a realization of kp. They dif- 
fer in how they compute these realizations and how they con- 
struct p(kp(kp)\y(ko'.kp)) from them, and, consequently, in 
their computational complexity. 

3.1. Monte Carlo Sampling 

To account for uncertainty in the prediction step, the simplest 
method is Monte Carlo sampling. Several realizations of the 


Algorithm 2 {k^}^ = MC(p(x(fc P ), 9(k P )\y(k 0 :k P )), 
p(X g ), p(X u ), p(X v ),N) 

1 : for i = 1 to N do 

2: (x w (fc P ), 9 {t) (kp)) ~ p{x(k P ), 9(kp)\y(k 0 :kp)) 

3: X ( g l) ~p(A e ) 

4: ®i P construct0(Ag l \ 0^\kp)) 

5: A„ ) ~p(A„) 

6: U« 4— constructU(A„*) 

7: Ai°~p(A„) 

8: Vj^ <— constructv(Aj, ! - > ) 

9: *£>«- P(xW(fcp) ) 0W,U«,V«) 

10: end for 


state, parameter trajectory, input trajectory, and process noise 
trajectory are sampled from their corresponding distributions. 
For each realization, kp is computed. The resulting set of kp 
values characterizes its distribution. 

Algorithm 2 shows the Monte Carlo prediction algorithm. 
The algorithm is given as input the joint state-parameter 
distribution, and the distributions of the Xg, A„, and A t , 
variables, along with the number of samples to take, N. 
For N times, the algorithm samples from the distributions, 
constructs the parameter, input, and process noise trajecto- 
ries, and calls the P function to compute kp. To construct 
the trajectories from the A variables, the construct© 
constructU and constructV functions must be pro- 
vided by the user, as they depend on the chosen surro- 
gate variables and how they are to be interpreted. Note 
that the construct© function requires an additional in- 
put, (kp), which is a sample from the estimator-computed 
joint parameter estimate at time kp. 

In any prediction algorithm, computational complexity is 
driven by two factors: the number of evaluations of P, and the 
length of time each sample takes to simulate to kp (Daigle & 
Goebel, 2010). Assuming a fair comparison for the second 
factor, we can compare algorithms mainly by the first factor. 
In the case of Monte Carlo sampling, the number of samples 
N is arbitrary and determines the efficiency. In most cases, a 
very large value of N is required in order to reproduce accu- 
rately the important characteristics of the kp distribution. 

3.2. Unscented Transform Sampling 

A more complex method that can improve the efficiency 
of prediction over the Monte Carlo method is to use the 
unscented transform (UT) to sample from the distribu- 
tions (Daigle, Saxena, & Goebel, 2012). We present here 
an extended and generalized version of the method originally 
presented in (Daigle, Saxena, & Goebel, 2012) in order to 
accommodate the A variable formulation. 

The UT takes a random variable a € R n " , with mean a and 
covariance P aa , that is related to a second random variable 
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b <E by some function b = g(a), and computes the 
mean b and covariance Pi,/, with high accuracy using a mini- 
mal set of deterministically selected weighted samples, called 
sigma points (Julier & Uhlmann, 1997). The number of sigma 
points is linear in the dimension of the random variable, and 
so the statistics (mean and covariance) of the transformed ran- 
dom variable can be computed more efficiently than Monte 
Carlo sampling. 

Here, A denotes the ith sigma point from a and w l denotes 
its weight. The sigma points are always chosen such that the 
mean and covariance match those of the original distribution, 
a and P aa . Each sigma point is passed through g to obtain 
new sigma points 13, i.e., 

B l = g (A 1 ) 


with mean and covariance calculated as 


b = ^ w 'B l 

i 

P b b = ^2w i (B i -b)(B i -h) T . 


In this paper, we use the symmetric unscented transform 
(SUT), in which 2 n a + 1 sigma points are symmetrically se- 
lected about the mean according to (Julier & Uhlmann, 2004): 


w 


A 1 


K 

(n a + k) ’ 

2 (n a + k) ’ 


i = 0 

* = !)•••) 2 n a 



where ^ \J (n a + K)P aa j refers to the zth column of the ma- 
trix square root of (n a + n) P aa , and n is a free parameter 
that can be used to tune higher order moments of the distribu- 
tion. When a is assumed to be Gaussian, (Julier & Uhlmann, 
1997) recommends setting n = 3 — n a . Note that with the 
UT, weights may be negative, and are not to be directly inter- 
preted as probabilities. 

For prediction, the G function serves as g in the above for- 
mulation, where a corresponds to the joint distribution of the 
state and A variables, and b corresponds to Ue- The pre- 
diction algorithm in this case is shown as Algorithm 3. The 
algorithm first uses the symmetric unscented transform to 
compute sigma points for the given probability distributions 
(treated together as a joint distribution), where each sigma 
point consists of a sample for the state-parameter vector and 
the A variables. For each sigma point, the parameter, input, 
and process noise trajectories are constructed and the P func- 


Algorithm 

3 {k^\wW}? =1 = UT (p(x(fc P ), 6{k P ) 

|y (k 0 :k P )),p(X e ), p{ A„), p{ A„)) 

1: 

N 4- 2 (r, 

ix+ng + n\ e + n\ u + n\ v ) + 1 

2: 

{x^(kp 

),0^{kp),\o,\u,\ v ,w (i) }? =1 <r- 


SUT((p(> 

i(kp),6{kp)\y(k 0 :kp)),p(\g),p(\ u ),p(\ v ))) 

3: 

for i — 1 

to N do 

4: 

©£ 4 

— construct0(Ag l \ 0^ l \kp)) 

5: 

u £>* 

- constructU(Au*) 

6: 

vj£ + 

- constructv(Aj, ! ' ) ) 

7: 


p(xW(fc P ),©W,U^,V«) 

8: 

end for 



tion is called to compute the corresponding Ice- The mean 
and variance for k can be computed from its sigma points. 

This prediction method will often require a smaller number of 
samples than with Monte Carlo sampling, since the number 
of sigma points grows only linearly with the problem dimen- 
sion. This is partly due to the fact that the UT method pro- 
vides only mean and covariance information, whereas addi- 
tional higher-order moments can be computed with the Monte 
Carlo method. Extended versions of the UT are also available 
that compute higher-order statistical moments (Julier, 1998). 

3.3. Inverse First-Order Reliability Method 

The Monte Carlo and UT approaches are sampling-based 
techniques to predict the uncertainty in the kp ■ Here we 
briefly explain an optimization-based method, the inverse 
first-order reliability method, for this purpose. The First- 
order Reliability Method (FORM) and the Inverse First-Order 
Reliability Method (inverse FORM) were originally devel- 
oped by structural engineers to evaluate the probability of 
failure of a given structure (Haidar & Mahadevan, 2000). In 
an earlier publication (Sankararaman et al., 2013), we have 
extended these two approaches for uncertainty quantification 
in the context of remaining useful life estimation, i.e., prop- 
agate the uncertainty in present estimates of states and pa- 
rameters, future loading, future process noise, and future pa- 
rameter values through P (defined earlier in Algorithm 1) to 
calculate the uncertainty in fc^. In the present paper, we use 
the inverse FORM methodology to calculate the entire prob- 
ability distribution of in terms of the cumulative distribu- 
tion function. Calculating the cumulative distribution func- 
tion (CDF) is equivalent to calculating the probability density 
function p(kE(kp)\y(ko-kp)), since the density function can 
easily be obtained by differentiating the cumulative distribu- 
tion function. 

For a generic random variable Z, the cumulative distribution 
function is a mapping from a realization 3 of the random vari- 
able to the set [0, 1], and is denoted by Fz(z). If Fz(z ) = p, 
then the probability that the random variable Z is less than 
a given value z is equal to //. In the context of prognostics, 
the goal is to compute the uncertainty in kp- Typically, the 
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Algorithm 4 fc_E(fcp) Pa(w) 

1 : 

[x(fcp) 

, 9(kp), \g, X u , A t , ] 4 — u 3 

2 : 

©ftp 5- 

- construct0(Ae, 9{kp)) 

3: 

U kp + 

- constructU(A„) 

4: 

Vfe p 4— constructv(A„) 

5: 

kE(kp 

) f- P(x(fcp), ©fc P , Ufcp, Vfcp) 


quantities x(fcp), 9(kp), A g, A„, and A„ are vectors, and 
now, consider a new vector which is the concatenation of all 
these vectors as ui = [x(fcp), 9(kp), A g, \ u , A„]. Based on 
the probability distributions of x(kp), 9(kp) A g, \ u , \ v , the 
joint probability density of u>, denoted as /o (w), can be easily 
calculated. Note that w is a realization of the random variable 
that is denoted by f 2 . 

In order to implement the inverse FORM, it is necessary to 
construct a function whose inputs are u = [x(fcp), 9(kp), A g, 
A„, A„] and the output is k p . This function similar to P, with 
one difference; while P takes realizations of entire trajecto- 
ries, i.e., ©fe^.Ufcp, and V kp as input arguments, the new 
function needs to consider realizations of the corresponding 
surrogate variables as input arguments. This new function, 
denoted by P^, is indicated in Algorithm 4. 

The inverse FORM approach is now explained using the func- 
tion k p = P x('jj). The reason for such vectorized represen- 
tation using io is not only to facilitate easy explanation of the 
FORM and inverse FORM algorithms but also demonstrate 
that these algorithms do not differentiate amongst state es- 
timate values, parameter values, future loading trajectories, 
and process noise trajectories but treat them similarly. 

Let the CDF of kp be denoted as Fp; e (fcp) = ??. Using 
kE = Pa(w), the FORM approach can be used to calculate 
the value of rj corresponding to a given value of kp. Con- 
versely, the inverse FORM approach can be used to calculate 
the value of /cp corresponding to a given value of //. By re- 
peating the FORM procedure for multiple values of fcp, or by 
repeating the inverse FORM procedure for multiple values of 
?/, the entire CDF F Kt (kp) can be calculated. In a practical 
scenario, it is not reasonable to know what values of fcp need 
to be selected to implement the FORM procedure, since the 
goal is actually to compute the uncertainty in kp. However, 
it is easier to select values of rj (say, 0.1, 0.2, 0.3 and so on 
until 0.9) which span the entire probability range and imple- 
ment the inverse FORM procedure for each of these 77 values. 
Therefore, we use the inverse FORM approach to quantify 
the uncertainty in fcp. The authors have explained the inverse 
FORM algorithm in detail in previous work (Sankararaman 
et ah, 2013); in this section, the overall approach is briefly 
summarized and the algorithm is provided. 

Both FORM and inverse FORM approaches linearize the 
curve represented by the equation fcp = P \(lu) and transform 
all the variables in il to standard normal variables (Gaussian 
distribution with zero mean and unit variance) using well- 


Algorithm 5 {k^\ »7^}£Lcl InverseFORM(/ n (w), Pa) 


1 : 

2 : 

3: 

4: 

5: 

6 : 

7: 

8 : 

9: 

10 : 

11 : 

12 : 

13: 

14: 

15: 

16: 

17: 

18: 

19: 

20 : 

21 : 


N 4 — Number of 77 values to consider. 

M 4— Number of elements in u 
{Example: N = 9, r / 1 - 1 = 0.1 x i, i = 1 to N. } 
for i = 1 to N (For every 77 value) do 
pW 4- — 4>— 1 ( t 7 ( * ) ) 

u>o Select initial guess for optimization 
Convergence = 0 
j = 0 {Iteration number} 
while Convergence 4 — 0 do 

4>j <— T(u>j) {Transformation to Std. Normal Space} 

<j>j [4>jk\ k = 1 to M] 

a.j <— [ajk; k = 1 to M\ where ajk = 7^77 


< t >3 + 1 



x P w 


uj.jp] T~ {</>j+ 1 ) {Transformation to Original Space} 
if ujj+i « ojj then 
Convergence 4 — 1 

end if 

3 <~j + 1 

end while 

k^E ^ ~ P A(Wj) 

end for 


defined, popular transformation functions. Thus, kE can be 
expressed as a linear sum of Gaussian variables, and therefore 
the probability distribution of kp can be computed easily. The 
key point in these algorithms is that the point of linearization 
is chosen to be the Most Probable Point (MPP), i.e., the point 
of maximum likelihood value. For example, in a Gaussian 
distribution, the MPP is at the mean. Since each uncertain 
quantity in il may have its own distribution, the MPP is com- 
puted in the standard normal space, where the origin has the 
highest likelihood of occurrence. However, the origin may 
not satisfy the equation kE = Pa(u>), and the point of liner- 
ization needs to lie on the curve represented by the equation 
kE = Pa(w). Therefore, the problem reduces to estimat- 
ing the minimum distance (measured from the origin, in the 
standard normal space) point on the curve represented by the 
equation kE = Pa(w). This is posed as a constrained min- 
imization problem, and solved using a well-known gradient- 
based optimization technique, as described in Algorithm 5. 
Once the minimum distance (denoted by /?) has been evalu- 
ated, then it can be proved that FK E {kE) = $(—/?), where 
$(.) represents the standard normal distribution function. Al- 
gorithm 5 explains the numerical implementation of the in- 
verse FORM approach. 

In the above algorithm, note that the user needs to specify 
functions T and T -1 for transforming original space to stan- 
dard normal space and from standard normal space to origi- 
nal space respectively. There are several types of transforma- 
tion available in the literature (Haidar & Mahadevan, 2000) 
and any valid transformation may be used. Further, note that 
the gradient a needs to be calculated in the standard normal 
space. This depends on ( 7 ) the gradient in the original space; 
and (if) the chosen transformation T. 

Using the algorithm, the values of kE corresponding to the 
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chosen values of i] are first obtained, and then an interpolation 
technique can be used to obtain the entire CDF. Also, if the 
goal is to quickly obtain bounds on Ue , then we may consider 
two 77 values in either tail of the probability distribution (say, 
for example, ?/ values of 0.1 and 0.9), and the corresponding 
probability bounds of can be obtained. 

4. Case Study 

As a case study, we perform battery prognostics on a plan- 
etary rover and present simulation-based results. The rover 
and its simulation are described in detail in (Balaban et al., 
2011). The rover battery system consists of two parallel sets 
of 12 batteries in series to provide around 48 V. We are inter- 
ested in predicting end-of-discharge (EOD), which is defined 
as the point when the voltage of a single battery drops below 
2.5 V. 

We consider two different scenarios for the rover. In both 
scenarios, the rover is provided a sequence of waypoints that 
must be visited. The rover travels through the waypoints in 
order until the batteries discharge. In the first scenario, the 
desired forward speed of the rover is the same when moving 
to each waypoint. In the second scenario, the desired forward 
speed is different depending on which waypoint is being ap- 
proached. Since the needed power draw from the batteries 
depends on speed, the first scenario resembles a constant- 
amplitude loading situation, and the second scenario resem- 
bles a variable-amplitide loading scenario. 

The battery prognostics architecture works as follows. The 
rover provides voltage measurements on all batteries, and the 
current going into the batteries. Because there are two sets of 
batteries in parallel, each battery sees only half the measured 
current. The measured current and voltage are fed into the 
battery prognoser. The battery prognoser uses an unscented 
Kalman filter (UKF) to perform state estimation (see (Julier 
& Uhlmann, 1997, 2004; Daigle, Saha, & Goebel, 2012) for 
details on the UKF). The state estimate is then fed into the 
predictor, which makes EOD predictions every 100 seconds. 


sp 1 L -S 



Figure 2. Battery equivalent circuit. 


tures the major nonlinear voltage drop due to surface overpo- 
tential, R s captures the so-called Ohmic drop, and R p mod- 
els the parasitic resistance that accounts for self-discharge. 
This empirical battery model is sufficient to capture the ma- 
jor dynamics of the battery while ignoring temperature effects 
and additional minor processes. The governing equations for 
the battery model are presented in continuous time below. 
The implementation of the proposed methodology considers 
a discrete-time version with a discrete time-step of 1 s. 

The state-of-charge, SOC, is computed as 

SOC = 1 - gm ° x ~ Qb , (5) 

^ max 

where (p, is the current charge in the battery (related to C)j, 
Qmax is the maximum possible charge, and C max is the max- 
imum possible capacity. The resistance related to surface 
overpotential is a nonlinear function of SOC : 

R'.sp ^spo exp ( R sp2 ( 1 - SOC)), (6) 

where R spo , R spi , and R SP2 are empirical parameters. The 
resistance, and, hence, the voltage drop, increases exponen- 
tially as SOC decreases. 


In the remainder of the section, we describe the details of the 
underlying battery model used by the rover simulation and the 
battery prognoser, and provide simulation-based experimen- 
tal results for different scenarios and comparing the different 
methods presented in Section 3. 

4.1. Battery Model 

The battery model is extended from that presented in (Daigle, 
Saxena, & Goebel, 2012). The model is based on an electrical 
circuit equivalent as shown in Fig. 2, following similar mod- 
eling approaches to (Chen & Rincon-Mora, 2006; Ceraolo, 
2000). The large capacitance Cb holds the charge (p, of the 
battery. The nonlinear Cb captures the open-circuit poten- 
tial and concentration overpotential. The R sp -C sp pair cap- 


Voltage drops across the individual circuit elements are given 
by 



(7) 

( 8 ) 

(9) 

( 10 ) 


where q sp is the charge associated with the capacitance C sp , 
and q s is the charge associated with C s . The voltage l), is 
also the open-circuit voltage of the battery, which is a nonlin- 
ear function of SOC (Chen & Rincon-Mora, 2006). This is 
captured by expressing Cb as a third-order polynomial func- 
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Table 1. Battery Model Parameters 


Parameter 

Value 

C’bo 

19.80 F 

c bl 

1745.00 F 

C b2 

-1.50 F 

c b 3 

-200.20 F 

Rs 

0.0067 Q 

Cs 

115.28 F 

Rp 

1 x 10 4 

C sp 

316.69 F 

Rsp 0 

0.0272 Q 

Rspi 

1.087 x 10“ 16 D 

R'SJ>2 

34.64 

Qmax 

3.11 x 10 4 C 

Cmax 

30807 C 


tion of SOC: 

C b = C bo + C bl SOC + C b2 SOC 2 + C b3 SOC 3 (11) 

The terminal voltage of the battery is 


V = V b -V sp - V s . 


( 12 ) 


Currents associated with the individual circuit elements are 
given by 


II 

(13) 

i b — ip P i , 

(14) 

. v sp 

*«p — p ) 

Rsp 

(15) 

. _ . v s 

IS t b jy • 

Rs 

(16) 

The charges are then governed by 


q b — i b t 

(17) 

Qsp Gp? 

(18) 

q s = i s - 

(19) 


In the case of the battery, the event E we are interested in 
predicting is EOD. Te is specified as V < Veod, where 
C E0D = 2.5 V. 

The parameter values of the battery model are given in Ta- 
ble 1. All voltages are measured in Volts, resistances are 
measured in Ohms, charges are measured in Coulombs, and 
capacitances are measured in Coulombs per Volt (or Farads). 
Note that C bo , C bl , C),, , and C b3 are simply fitting parameters 
in Eq. 1 1 and do not have physical meaning. 

For the battery model, x = [q b q sp q s ], 6 = 0 (i.e., all pa- 
rameters are assumed constant and no parameters will be es- 
timated online), and y = [ V ]. We consider power P to be 
the input to the battery, so i = P/V, i.e, u = [/’]. Here, we 
choose power as the input, rather than current as in previous 
battery prognostics works, because it is simpler to describe 


battery load in terms of power. For the same power demands 
from the rover onto the battery, current will increase as bat- 
tery voltage decreases; it is necessary to capture this current- 
voltage relationship in order to use current as input. There- 
fore, it is much easier to predict future power usage than to 
predict future current draw, and hence, power is used as input. 

4.2. Future Input Characterization 

In the experiments presented in this section, we will consider 
only uncertainty in the state and in the future inputs (meth- 
ods for dealing with process noise are described in (Daigle, 
Saxena, & Goebel, 2012; Sankararaman & Goebel, 2013)). 
Therefore we need to define p(X u ) and the constructU 
function. We explore three methods that differ in complexity 
and the amount of system knowledge used. 

The future input trajectory to a battery model depends on how 
the rover will be used. When moving from one waypoint to 
the next, the rover must turn towards the next waypoint while 
maintaining its foward speed (that is how the locomotion con- 
troller is designed to work). For the same forward speed, 
turning actually requires more power than going straight, be- 
cause the rover must also move against the ground torques 
produced while turning in addition to the opposing forces pro- 
duced when moving forwards. Further, because of the turn- 
ing, the actual distance traveled between two waypoints is 
greater than the straight-line distance between them because 
the rover actually takes a curved path. 

To correctly account for all these details, a system-level ap- 
proach is required (Daigle, Bregon, & Roychoudhury, 2012). 
In this case, the whole rover and its locomotion controller 
would be considered as the system under prognosis. Thus, the 
whole rover would be simulated moving through the different 
waypoints, and this would define very precisely (depending 
on model fidelity) the power drawn from the batteries as a 
function of time. However, such an approach is more com- 
putationally expensive than considering only a single battery 
model. 

A simpler approach is to assume that, in the current opera- 
tion of the system, the future inputs to the system will look 
like the past inputs. That is, we can assume that the future 
power requirements will be constant with the mean and vari- 
ance defined by the past power requirements over some time 
window. If the window size is large enough, then the differ- 
ences in power that arise may be represented well enough in 
the statistics of the past behavior. Although simple, such an 
approach may work well in some circumstances. 

As a middle ground, we can incorporate some system-level 
knowledge into predicting the future power requirements 
without resorting to a system-level simulation. We can do this 
by computing the mean and variance of the power draw and 
distance traveled between waypoint pairs for each forward 
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speed setting. Then, knowing the current rover location and 
the remaining waypoints, given a realization of each distance 
and power variable for the remaining waypoints and the de- 
sired forward speed when heading to each waypoint, we can 
compute the future power as a function of time based on the 
expected path the rover is going to take. This approach uses 
system knowledge, i.e., knowledge beyond the battery model, 
in order to compute useful predictions of the future inputs to 
the battery, and therefore makes useful EOD predictions for 
the battery, but without resorting to a system-level simulation. 

4.3. Results 

In the simulation experiments considered in this section, all 
parameters are considered known exactly and no process 
noise is added. The two potential sources of uncertainty are 
related to the state estimate obtained by the UKF and the fu- 
ture input assumptions. Predictions are made every 100 s un- 
til EOD, and the accuracy and precision metrics are averaged 
over all these predictions. We use the relative accuracy (RA) 
metric as defined in (Saxena, Celaya, Saha, Saha, & Goebel, 
2010) as a measure of accuracy and relative standard devia- 
tion (RSD) as a measure of spread. In the following plots, the 
* superscript indicates the ground truth values. 

4.3.1. Constant-Loading Scenario 

We consider first the scenario where the rover must move be- 
tween equidistant waypoints at the same forward speed, re- 
sembling a constant-loading situation. Let us first assume that 
the future inputs (the battery power) are known exactly. There 
are 3 states in the battery model, so 7 sigma points are used 
by the UKF, and these are directly simulated forward to com- 
pute the Ue distribution using the sigma point weights. In 
this case, since uncertainty is limited only to that in the state 
estimate, predictions are both very accurate and precise, with 
RA averaging to 99.65% and RSD to 0.64%. 

Now assume that the past power requirements are statistically 
representative of the future power requirements. Here, we 
consider Monte Carlo sampling with 3500 samples. We con- 
sider window sizes of 100 s, 500 s, and unlimited size. The 
results are shown in Fig. 3. In all three cases, the uncertainty 
starts initially very large, because the window is not large 
enough to accurately capture the statistics of the power usage. 
With a small window size (Fig. 3a), the statistics of the power 
usage averaged over the window fluctuate. The variance is 
larger when both turns and forward movements appear in the 
same window, and smaller when only forward movements are 
in the window. With a larger window size the variance will 
average to a larger value that accounts for both turns and for- 
ward movements, as seen in Figs. 3b and 3c. Since the past 
power usage turns out to be a good indicator of future power 
usage in this scenario, the results are fairly accurate, with RAs 
of 98.14%, 98.47%, and 97.53% for 100 s, 500 s, and un- 
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(a) Window size of 100 s. 
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Figure 3. A kE predictions using Monte Carlo sampling. 


limited window sizes, respectively. Corresponding RSDs are 
3.83%, 5.57%, and 7.61%. The spread increases as window 
size increases since more variation is accounted for in larger 
windows. 

Using a window size of 500 s, the results using the UT are 
shown in Fig. 4. Here, the results are comparable to using 
Monte Carlo sampling with the same window size, with an 
average RA of 98.69% and RSD of 5.24%. The UT method, 
however, needs only 9 total samples, with there being only 
3 states and one surrogate input variable to consider. This 
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Figure 4. A predictions with a window size of 500 seconds 
using UT sampling. 
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Figure 5. Afce predictions with a window size of 500 seconds 
using inverse FORM. 

results in a substantial decrease in computational cost com- 
pared to the Monte Carlo approach. Fig. 5 shows the results 
using inverse FORM. The results are similar to using both 
Monte Carlo and UT sampling, and RA averages to 98.46% 
and RSD to 3.43%. 

Using knowledge about the future waypoints to be visited, we 
can improve over using a window of past data to determine 
future inputs to the system. Fig. 6 shows the improved future 
input characterization method using Monte Carlo sampling 
with 3500 samples. The plots look the same for UT sam- 
pling and inverse FORM. The accuracy is comparable to the 
previous approach, with an average RA of 98.29% for Monte 
Carlo sampling, 98.32% for UT sampling, and 98.30% for in- 
verse FORM. RSDs, however, are lower now since the future 
inputs are known with more precision than could be derived 
from a window of past samples. RSD averages to 1.61% for 
both Monte Carlo and UT sampling and 4.67% for inverse 
FORM. 

4.3.2. Variable-Loading Scenario 

We now consider the second scenario that uses the same way- 
points as the previous scenario, but the rover is commanded 
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Figure 6. A Ue predictions with improved future input char- 
acterization using Monte Carlo sampling. 


to go different speeds depending on which waypoint is being 
headed towards, resembling a variable-loading situation. As- 
suming the future inputs are known exactly, the average RA is 
99.50% and the average RSD is 0.70%. The only uncertainty 
is in the state estimate. 

Because the speed of the rover will change with each new 
waypoint, it is no longer correct to assume that past power re- 
quirements are representative of future power requirements. 
Fig. 7a shows the results when we incorrectly make this as- 
sumption for Monte Carlo sampling with 3500 samples, us- 
ing a window size of 500 s. Clearly, the predictions are very 
inaccurate. RA averages to 90.53% and RSD to 13.74%. Us- 
ing UT sampling we find similar results, with RA averaging 
to 90.12% and RSD to 13.54%. Using inverse FORM, RA 
averages to 90.59% and RSD to 9.24%. When the average 
speed of the rover in the future is greater than what is as- 
sumed based on the window of past samples, then A &e is 
overestimated. When the average speed is less than what is 
assumed, A Ue is underestimated. Because the average speed 
over the window changes based on the previous waypoints 
within that window, the predictions fluctuate above and be- 
low ground truth. If the window size is increased, such that 
it accounts for all possible speed settings, then accuracy can 
be improved because the assumed average speed based on 
past samples will match the future average speed, however, 
spread will also increase since multiple speeds are consid- 
ered in the window. Predictions with an unlimited window 
size are shown in Fig. 7b. The predictions initially fluctuate 
as the window begins to fill up, but by 2000 s the predictions 
have smoothed out. The spread is clearly larger than with the 
smaller window size, but predictions are more accurate once 
the window contains all potential future speeds. In this case, 
RA improves to 96.62% but RSD increases to 17.47%. 

Predictions can be improved by using system knowledge to 
help characterize the future inputs. In this case, the future 
power as a function of time is computed based on the rover’s 
current location, the remaining waypoints, and the desired 
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Figure 7. A Ice predictions using Monte Carlo sampling. 


speeds when heading to each waypoint, and measured statis- 
tics on average power between waypoints and average dis- 
tance traveled (to account for turns). Fig. 8 shows the results 
using Monte Carlo sampling. The plots for UT sampling and 
inverse FORM look the same. All three methods are now 
clearly very accurate and precise. RA averages to 98.85% 
and RSD to 1.95% for Monte Carlo sampling, 98.82% and 
1.91% for UT sampling, and 97.77% and 2.72% for inverse 
FORM. Unlike when using a window, here knowledge of the 
future waypoints and desired speeds allows accurate predic- 
tions to be made from the start of the scenario, and with very 
little spread. That is, the future inputs are well-known and so 
predictions are very close to the optimal. 

In the case above, for UT sampling, there are 3 states to con- 
sider and at most 100 surrogate input variables. For the in- 
put variables, there are two variables associated with each re- 
maining waypoint, one for the power that will be consumed 
heading towards the waypoint and one for the distance to 
travel to a waypoint from the previous waypoint (due to turn- 
ing while moving forward, the distance is more than the linear 
distance between the waypoints). Since there are 50 way- 
points, there are 100 random variables needed. This yields 
2(103) + 1 = 207 samples, which is relatively small com- 
pared to what Monte Carlo sampling would require to achieve 
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Figure 8. A predictions with improved future input char- 
acterization using Monte Carlo sampling. 


the same performance. As the rover visits waypoints, the 
number is reduced, so 207 samples is only the maximum. 
When only a single waypoint is left, only 2(5) + 1 = 11 
sigma points are needed. 

4.4. Discussion 

The two scenarios demonstrate the importance of the future 
input characterization problem. Even though the rover is a 
complex system, in the first scenario, the simple assumption 
that the future inputs will look like the past inputs was suf- 
ficient for accurate and precise predictions. The additional 
power required by turns was captured using the statistics 
of a window of past samples, so that they did not have to 
be explicitly accounted for and assuming constant future in- 
puts was sufficient. Using system-level information about the 
waypoints the rover would visit improved significantly on the 
uncertainty associated with the future inputs but did not sig- 
nificantly impact accuracy. 

In the second scenario, the assumption that the future inputs 
look like the past inputs did not provide as accurate or pre- 
cise results as with the first scenario. Performance could 
have been potentially worse if the rover did not cycle through 
different speed settings and instead, for example, always in- 
creased the speed for the next waypoint. In this scenario, truly 
accurate predictions could be made only when system-level 
knowledge was utilized to predict the battery inputs. This ap- 
proach still made some simplifying assumptions allowing a 
component-level prognostics approach to still be used. 

Given a particular method for future input characterization, 
Monte Carlo sampling, UT sampling, and inverse FORM all 
had comparable accuracy and precision. With Monte Carlo 
sampling, a large number of samples were used and with 
smaller numbers of samples, performance decreases, so the 
number of samples required depends on the prognostics per- 
formance requirements. In this sense Monte Carlo sampling 
has an advantage because its computational complexity can 
be tuned. In addition it is the relatively simplest approach 


11 



Annual Conference of the Prognostics and Health Management Society 2013 


to implement of the three methods. A disadvantage is that 
it is a stochastic algorithm, which can be problematic for 
verification and certification procedures (Daigle, Saxena, & 
Goebel, 2012). UT sampling, in contrast, is a determinstic 
algorithm, and it selects only the minimal number of sam- 
ples and this number grows only linearly with the number 
of random variables. A disavdvantage though is that it com- 
putes only mean and variance of the predictions, which may 
not be enough information in some cases. Inverse FORM, 
on the other hand, is not only deterministic but also allows 
control of both computational complexity and accuracy by 
selecting desired CDF values and computing the correspond- 
ing percentile values of Ate- The probability distribution of 
Jze can be reconstructed from that information and any de- 
sired statistical moments may be calculated. For each inverse 
CDF calculation, three to four iterations are usually required 
for optimization convergence. If the number of random vari- 
ables is denoted by n (length of vector u), each iteration of 
Inverse-FORM requires n + 1 sample evaluations of G, where 
one evaluation is required for computing kpj and n evalua- 
tions for computing the gradient vector of ks- Therefore, if it 
is desired to repeat inverse FORM for k different CDF values, 
then k x 4 x (n + 1) evaluations of P are required. Thus, the 
computational complexity linearly increases with the number 
of random variables and results in increased information re- 
garding uncertainty. 

5. Conclusions 

In this paper, we provided a general formulation of the prog- 
nostics problem and its uncertainty. Given descriptions of 
the sources of uncertainty, i.e., state uncertainty, parameter 
uncertainty, future input uncertainty, and process noise, we 
provided an algorithmic framework for incorporating this un- 
certainty into the predictions. With the novel concept of sur- 
rogate variables, we presented three methods for propagating 
the uncertainty: Monte Carlo sampling, unscented transform 
sampling, and the inverse first-order reliability method. Using 
battery prognostics on a planetary rover as a case study, we 
proposed two future input characterization methods and com- 
pared the performance of the different prediction algorithms 
for these methods for different scenarios in simulation. All 
approaches had similar performance, yet each offer different 
advantages and disadvantages that suggest when one would 
be preferred over another. 

In future work, we will further investigate these ideas on other 
systems, and further develop the uncertainty quantification 
framework. While the proposed methods are promising for 
estimating uncertainty in prognostics, their applicability to 
multi-modal probability distributions, particular in the con- 
text of remaining useful life estimation, needs to be investi- 
gated. Further, we will also focus on model uncertainty quan- 
tification and develop methods for estimating model errors 
and model parameter uncertainty separately, instead of sim- 


ply using lumped process noise terms. 
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